Tunneling of Bound Systems at Finite Energies: Complex Paths Through Potential 

Barriers. 
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We adapt the semiclassical technique, as used in the context of instanton transitions in quantum field 
theory, to the description of tunneling transmissions at finite energies through potential barriers by 
complex quantum mechanical systems. Even for systems initially in their ground state, not generally 
describable in semiclassical terms, the transmission probability has a semiclassical (exponential) 
form. The calculation of the tunneling exponent uses analytic continuation of degrees of freedom 
into a complex phase space as well as analytic continuation of the classical equations of motion into 
the complex time plane. We test this semiclassical technique by comparing its results with those of 
a computational investigation of the full quantum mechanical system, finding excellent agreement. 
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1. Tunneling phenomena are inherent in numerous 
quantum systems, from atoms to condensed matter to 
quantum field theory. Even in systems with a small 
parameter — coupling constant — a quantitative descrip- 
tion of tunneling is possible only in a limited number 
of cases. Perhaps the best known example is the WKB 
approximation familiar from one-dimensional wave me- 
chanics; similar techniques, such as the "most probable 
escape path" and instanton methods [0J^], are used to 
study tunneling from the bottom of potential wells. In 
the latter cases the calculation of the tunneling proba- 
bility may be reduced to the solution of classical equa- 
tions of motion for real generalized coordinates in imag- 
inary ("Euclidean") time, supplemented by the analysis 
of small fluctuations about this classical Euclidean tra- 
jectory. However, these methods often fail in describing 
tunneling of systems with more than one degrees of free- 
dom at finite energies. 

It has been suggested recently [||^], in the context of 
instanton transitions in quantum field theory, that semi- 
classical techniques may be used for calculating the expo- 
nential suppression factors in a class of processes where 
multi-dimensional systems tunnel at finite energies. The 
proposal involves a double analytic continuation: the 
degrees of freedom are continued into a complex phase 
space, and the equations of motion are solved along a 
contour in complex time. The tunneling exponent is de- 
termined by an appropriate solution of the calssical, al- 
beit complexified, equations of motion. Computation by 
numerical methods is then feasible even for systems with 
a large number of degrees of freedom, as has already been 
demonstrated in a field theoretic model M. A problem 
with the formalism of Refs. |||| is that its derivation 
from first principles is still lacking, although its plausi- 
bility has been supported by perturbative calculations 



about an instanton 

The purpose of this paper is two fold. First, we adapt 
the technique of Refs. to tunneling of quantum me- 
chanical bound systems through high and wide potential 
barriers. As an example, we consider a system of two 
degrees of freedom with linear binding force. We find 
that if the bound system is initially in a highly excited 
state, the tunneling exponent is indeed calculable in a 
semiclassical way. This result is hardly surprising, as the 
initial state itself can be described in semiclassical terms. 
We formulate the complexified classical boundary value 
problem relevant to the calculation of the exponent in 
this case. 

Second, the real strength of this formalism is that it 
also enables one to treat barrier penetration when the 
bound system is initially in a low lying state, e.g. the 
ground state. This is far from obvious, as this initial 
state cannot be described semiclassically. Nevertheless, 
we argue that in this case the tunneling exponent can 
be obtained by an appropriate limiting procedure. The 
resulting technique is less-well justified, so we have cho- 
sen to test it by direct computation of the transmission 
probability in the full quantum theory. We briefly de- 
scribe the numerical methods involved, and present the 
results of both the full quantum mechanical and semiclas- 
sical analyses. We find good agreement between the two, 
confirming the validity of the semiclassical approach. 

2. To be specific, let us consider a quantum mechanical 
system of two particles of equal mass m = 1/2 moving in 
one dimension. Let these particles be bound by the har- 
monic potential (uj 2 /8)(xi — X2) 2 , and one of these parti- 
cles be repelled from the origin by a positive semidefinitc 
potential V{x\) that vanishes as x\ — ► ±00 (we could 
of course allow V to depend on X2 as well, provided it 
couples to the internal degree of freedom). We take this 
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potential to have the form V(xi) = g~ 2 U(gxx), where g 
is a small constant. We set % = 1, so the classical limit 
corresponds to g — > 0. In what follows we present the re- 
sults of numerical calculations for u> = 1/2 and gaussian 
potential, t/(x) = cxp(— x 2 /2), although the treatment of 
other potentials would be similar. In terms of the center- 
of-mass and relative coordinates, X = (x\ + X2)/2 and 
y = [xi — x%)/2, the Lagrangian reads 

L = \x 2 + X -y 2 - ^y 2 - ±U[g(X + y)} (1) 

Far from the origin (\X\ — > oo), the center-of-mass and 
internal degrees of freedom decouple and the system can 
be characterized by its center-of-mass momentum P and 
oscillator excitation number n, or, equivalently, by n and 
the total energy E = P 2 /2+w(n+l/2). We wish to calcu- 
late the probablity T n (E) for transmission of the system 
through the barrier V. Of particular interest is Tq{E), 
the transmission probability of this system initially in its 
oscillator ground state. 

It is convenient to introduce rescaled total energy and 
occupation number e = g 2 E and v — g 2 n. With our 
choice of U, the top of the barrier corresponds to a poten- 
tial energy e = 1. For e < 1 transmission is possible only 
via tunneling. For e just above 1 classical over-barrier 
transitions are possible for very special initial states. In- 
deed, there exists an unstable, static classical solution 
with both particles stationary at the top of the barrier, 
x\ = X2 = 0, so that e = 1. If one perturbs this solution 
by giving an arbitrarily small, common positive velocity 
to both particles, they will move toward X = oo. The re- 
versed evolution takes the system to X = — oo, with the 
classical oscillator characterized by a certain excitation 
energy eg sc = cj^o (and a certain phase of the classical 
oscillator). The combined evolution is the classical tran- 
sition over the barrier from this particular asymptotic 
state. By solving the (real time) classical equations of 
motion numerically, we found that vq ~ 0.9 for uj = 1/2. 

The classical evolution of the system initially in the 
classical oscillator ground state {x\ = X2) leads to the 
excitation of the oscillator as it approaches the barrier. 
Classical transition over the barrier occurs in this case 
only if the total energy exceeds some critical value. In 
our example we found numerically e cr u =1.8. 

If e and v are such that classical transitions over the 
barrier are not possible, the system has to tunnel. We will 
shortly see that at e and v fixed, and g — > (i.e., at large 
total energy and initial occupation number, E, n oc g~ 2 ) 
the transmission probability has the semiclassical form 

T n {E) = C{e,v)e~^ F{t ' u) (2) 

The case of the initial oscillator ground state is more sub- 
tle. In analogy to Refs. |||| we suggest that the trans- 
mission probability at n — has the form 



T (E) = C (e)e-^ Fo{e) (3) 

and that the exponent is obtained by taking the limit 

F (e)=]hn v ^ F(e,u) (4) 

One of the main purposes of this paper is to check this 
limiting procedure by comparison with a fully quantum 
mechanical calculation. 

3. To see that Eq.(j|) is indeed valid, and to obtain 
the procedure for calculating the exponent F(e, v), let us 
consider the transmission amplitude A(Xf,yf,P,n) = 
(Xf,y f \exp[-iH(t f -ti)]\P,ri), where X f (> 0) and y f 
are the coordinates at time t / , and we eventually take the 
limit (tf — ti) — ► 00. This amplitude may be written as 
a convolution of the evolution operator in the coordinate 
basis and the wave function of the initial state. The for- 
mer is given by the path integral (Xf,yf \ exp[— iH{tf — 
ti)]\Xi,yi) = J[dX][dy] cxpiS where the integration runs 
over paths satisfying (X,y)(U) = pQ,^), (X,y)(tf) = 
(Xf,yf). For an initial state with P oc g^ 1 , n oc g~ 2 , 
the initial wave function is semiclassical and has the ex- 
ponential form. In the case of a harmonic binding poten- 
tial, this follows from the integral representation in the 
coherent state formalism: 

(X uyt \P,n) = ^==- A^ e -««-^e-*-*-W4VS»« 
V27T J 2m Vn! 

(One may replace z n j\fn\ by exp(n log z/y/n + n/2) 
at large n.) By introducing the rescaled integration 
variables X — > gX, y — > gy, etc., we observe that 
A(Xf,yf,P,n) is given by an integral of an exponential 
of the form exp(— g~ 2 T) where T depends only on the 
rescaled integration variables, v and e, and does depend 
explicitly on g 2 . This allows for a semiclassical analysis: 
we find stationary points of T and evaluate the integrals 
using a stationary phase approximation. We outline the 
main steps in the derivation of the stationary point equa- 
tions. 

Variation of V with respect to X(t) and y(t) for ti < 
t < tf leads to the conventional classical equations of 
motion. When classical transitions are forbidden, there 
will be no real solutions satisfying the boundary condi- 
tions. Nevertheless there will be solutions with complex 
values of the integration variables. When performing the 
analytic continuation we will, in general, encounter sin- 
gularities. To deal with this problem, we note that the 
time contour, originally the real axis, can be distorted 
into the complex plane, keeping the end points ti, tf 
fixed. This deformation of the time contour allows us to 
avoid these singularities. Thus, our strategy is to search 
for complex solutions of the classical equations of motion 
along a contour ABCDE in the complex time plane, as 
shown in Fig. 1. 
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There are further stationary point equations coming 
from variation of T with respect to the integration vari- 
ables at the end point i,-. It is convenient to formu- 
late these equations along part B of the contour, where 
t = iT/2 + t', t' = real — * — oo (this is possible because 
the equations of motion decouple in the asymptotic past). 
Instead of e and v we introduce new real parameters T 
and 8; T enters the problem through the shape of the 
contour. The general complex solution at large negative 
t' is X(t') = X + pt', y{t') = ue- wt ' + ve lut ' where 
Xq, p, u and v are complex parameters. The stationary 
point equations at the initial time lead to the following 
boundary conditions: (i) X(t') is real (i.e. p is real and 
T may be chosen so that Xq is also real) , (ii) the positive 
and negative frequency parts of y(t') are related to each 
other by v = u*e e . 

More boundary conditions appear when one evalu- 
ates the total transmission probability, i.e. integrates 
\A(Xf,yf-,P,n)\ 2 over Xf and y/, again in a gaussian 
approximation. These conditions involve the final time 
and simply require that (hi) X(t) and y(t) are real on 
the DE part of the contour. 



Im t 



® 



^> Re t 



Fig 1. Complex time contour used to find the stationary 
point solutions. 

At given T and 9 these three boundary conditions are 
sufficient to specify the complex solution of the classical 
equations of motion on the contour BCDE (up to time 
translations along the real axis). Given this solution, the 
exponent for the transmission probability (|2|) is the value 
of 2Rer at the stationary point. Explicitly, we find 



F(e, v) = 2ImSo - eT - v6 



where 

So = — 



dt 



BCDE 
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-Xd(X + -yd(y + —y' + U(X + y) 



is the (rescaled) classical action for the complex solution 
of the above boundary value problem. The total energy 
and excitation number are related to T and 6 by 



d(2Im5o) 
df 



0(2ImS o ) 



i.e. the pairs (e, v) and (T,ff) are Legendre-conjugate. 

We have solved the equations of motion numerically 
along the contour BCDE subject to the boundary condi- 
tions (i)-(iii). In particular, we have evaluated the limit 
(Q). The result of this semiclassical calculation is shown 
in Fig. 3. 

4. To check this semiclassical procedure, we have per- 
formed a numerical analysis of the full quantum system 
defined by (|l]). This is conveniently done in a basis of 
center-of-mass coordinate X eigenstates and oscillator 
excitation number n. In this basis the state is represented 
by a multi-component wave function ip n (X) = (AT, rz | M/) , 
and the time-independent Schrodinger equation reads 
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UJlj} n (X) 



dX 2 1 V" ' 2, 

+ ^Vnn'{X)lpn'{X) = Ehj) n {X) (5) 
n' 

where V nn '(X) = (n\V(X + y)\n'). Our choice of a gaus- 
sian potential V enables us to calculate V nn '(X) by a 
numerical iteration procedure. Equation (^|) is supple- 
mented with the standard boundary conditions: (a) the 
incoming wave (X — + — oo) is in a state of given center- 
of-mass momentum P and excitation number n; (b) only 
outgoing waves exist at X — * +oo. 

To solve the system (||) numerically, we introduce a lat- 
tice with equal spacing, X% — ka, and discretize eq. (|5|) 
using the Numerov-Cowling algorithm (which reduces 
the discretization error to 0(a 6 )). We also truncate the 
system to a finite number of oscillator modes n < Nq. In 
order to insure good accuracy of the solution, we have 
chosen the number of lattice sites 2Nx and the cutoff A*o 
as large as 2N X = 2-4096, N Q = 400. This corresponds to 
over 3 million coupled complex equations. To deal with 
them, we take advantage of the special form of Eq. (jjj). 
Indeed, by inverting a set of (A^o + 1 ) x (Nq + 1) matrices, 
which is computationally feasible, Eq. (g) can be recast in 
the form ip n (X k ) = En'[^^'(^-i) + Rkipn'(X k +i)}. 
The elimination of ip n at definite X k leads to a system 
of similar form for the remaining variables (with suit- 
ably redefined L and R), again after (No + 1) x (Nq + 1) 
matrix algebra and matrix inversion. In this way we pro- 
gressively eliminate variables at intermediate values of 
Xk and ultimately obtain a system that linearly relates 
ip n at the end points X — —Nxa and X = +Nxa. With 
a discretized version of the boundary conditions (a) and 
(b), this final system is straightforward to solve. The 
transmission probability is then determined by \ipn\ 2 at 
the end point X = N x a. 

We performed a series of checks of this numerical pro- 
cedure to insure that our calculations are sufficiently pre- 
cise and that the results are close to the continuum limit. 
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Fig 2. Logarithm oj 'the transmission probability vs. l/g 2 . 

We present in Figs. 2 and 3 the results of the full quan- 
tum mechanical computation of the transmission proba- 
bility for the system initially in its oscillator ground state. 
The potential V is gaussian, and oj = 1/2. Figure 2 shows 
that the transmission probability Tq(E) indeed has the 
functional form (||): at fixed e = g 2 E, the logarithm of 
Tq is very well fit by a linear function of g~ 2 . We use this 
fit to obtain the exponent F (e). Both the full quantum 
mechanical results for Fo(e) and the semiclassical results 
(the latter obtained by implementing the limiting pro- 
cedure (|])) are shown in Fig. 3. Clearly there is good 
agreement between the two. (The slight discontinuities 
in the quantum mechanical results are an artifact of the 
energy dependence of the g 2 range from which we can 
extract Fq. They provide an indication of the errors due 
to higher order effects.) We conclude that the validity 
of the semiclassical approach is confirmed by the direct 
quantum mechanical computation. 




Fig 3. Quantum mechanical and semiclassical results. 

5. Full quantum mechanical computations (analytic or 
numerical) of barrier penetration probabilities are rarely 
possible. Even for our simplified system, values of g 
smaller than 0.1 are difficult to study as one has to deal 
with very small transmission coefficients. On the other 



hand, limitations of the semiclassical computations are 
far less severe. The generalization of the semiclassical 
approach to quantum-mechanical systems with harmonic 
binding of more than two particles in more than one 
space dimension is straightforward, and we also expect 
that other binding potentials may be treated in a sim- 
ilar way provided their semiclassical wave functions are 
known. Indeed, in all such cases the transmission ampli- 
tudes with highly excited initial states will be given by 
(path) integrals of exponential functions, and the tunnel- 
ing exponents will be determined by appropriate station- 
ary points. The latter will be complex solutions to clas- 
sical field equations on contours in complex time, with 
boundary conditions depending on the binding potential. 
A limit analogous to Eq. (^) will then determine the tun- 
neling exponent for incoming systems in low lying bound 
states. 

The semiclassical calculability of pre-exponential fac- 
tors is less clear. While it is plausible that these fac- 
tors are given by functional determinants about complex 
classical solutions for highly excited incoming states (fi- 
nite v in our model), we do not expect that a limiting 
property similar to Eq. (Ej) will continue to hold for the 
pre-exponents. The calculation of such pre-exponential 
factors for low lying states remains an interesting open 
problem. 
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